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PROPER ORTHOGONAL DECOMPOSITION IN OPTIMAL CONTROL OF FLUIDS 


S.S. RAVINDRAN' 

Abstract. In this article, we present a reduced order modeling approach suitable for active 
control of fluid dynamical systems based on proper orthogonal decomposition (POD). The rationale 
behind the reduced order modeling is that numerical simulation of Navier-Stokes equations is still 
too costly for the purpose of optimization and control of unsteady flows. We examine the possi- 
bility of obtaining reduced order models that reduce computational complexity associated with the 
Navier-Stokes equations while capturing the essential dynamics by using the POD. The POD allows 
extraction of certain optimal set of basis functions, perhaps few, from a computational or experimen- 
tal data-base through an eigenvalue analysis. The solution is then obtained as a linear combination 
of these optimal set of basis functions by means of Galerkin projection. This makes it attractive for 
optimal control and estimation of systems governed by partial differential equations. We here use it 
in active control of fluid flows governed by the Navier-Stokes equations. We show that the resulting 
reduced order model can be very efficient for the computations of optimization and control problems 
in unsteady flows. Finally, implementational issues and numerical experiments are presented for 
simulations and optimal control of fluid flow through channels. 

Key words. POD, reduced order model, flow control, optimal control, Galerkin methods. 

AMS subject classifications. 93B40, 49M05, 76D05, 49K20, 65M60, 76D15 

1. Introduction. The invention of Micro Electro Mechanical Systems and other 
fast micro-devices has generated substantial interest in active control of fluid dy- 
namical systems for the design of advanced fluid dynamical technology. There are 
a large number of articles devoted to this actively growing field. For example, in 
[7, 9, 8, 10, 26, 4] various optimal control problems in viscous incompressible flows 
were discussed. In [22, 17, 21, 16, 25] experimental efforts were reported. However, ef- 
ficient computational methodologies for use in on-line, real-time computation for PDE 
based control design has seen little progress. In this article we discuss a reduced order 
method for PDE based control using the proper orthogonal decomposition (POD). 

The solution of complex fluid dynamic equations using the available finite element, 
finite volume, finite difference or spectral methods is, in general, not feasible for real- 
time estimation and control. There are methods that would yield small degree of 
freedom models for the purpose of control of partial differential equations. However, 
they do not adequately represent the physics of the system and may be very sensitive 
to operating conditions as they are based on input/output data of a given system. 

We here examine the possibility of obtaining reduced order models that reduces 
computational complexity associated with the Navier-Stokes equations while captur- 
ing the essential dynamics by using the proper orthogonal decomposition (POD). The 
proper orthogonal decomposition is a model reduction technique for complex nonlin- 
ear problems. It was first, proposed by Karhunen [11] and Loeve [14] independently 
and sometimes called Karhunen-Loeve (K-L) expansion. Subsequently it has been 
applied in various applications. In [15] the method was first called POD and there it 
was used to study turbulent flows. In [23] another important progress was made and 
the method of “snapshots” was incorporated into the POD framework which will be 
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described in the sequel. Other applications in turbulent flow simulations are given in 
[1, 24, 2, 18, 3, 20] and [12], 

When discretizing nonlinear partial differential equations using finite volume, fi- 
nite difference, finite element or spectral methods, one uses basis functions that have 
very little connection to the problem or to the underlying partial differential equations. 
In some spectral methods Legendre polynomials are used, in finite element methods 
piecewise polynomials are used and in finite difference methods grid functions are 
used. However, POD uses basis functions that are generated from the numerical 
solutions of the system or from the experimental measurements. 

The essential idea is to generate optimal basis functions for Galerkin representa- 
tions of PDEs. In other words, given an ensemble S = {ULljAj consisting of N data 
vectors of length N s , the K-L theory yields that we can find an orthonormal coordi- 
nate system {V^)}^ that is optimal in the sense that the variance of the dataset in 
the coordinate directions becomes maximal. Thus, when the Navier-St.okes equations 
are projected onto this optimal base using a Galerkin projection, one obtains a re- 
duced order model. The beauty of the POD is that it is a nonlinear model reduction 
approach and its mathematical theory is based on the spectral theory of compact, self 
adjoint operators. 

Our goal here is to discuss a computational approach based on reduced order 
models resulting from the application of POD for the active control problems arising 
in nonlinear fluid dynamic systems. 

As a test problem we take a two-dimensional flow through backward facing step 
channel. This flow configuration is considered as a typical unsteady separated flow. 
For high Reynolds’ numbers, flow separates and recirculation appears. We will for- 
mulate a recirculation control problem in this configuration with the control action 
achieved through the surface movement/blowing of mass on a part of the boundary. 

The plan of the paper is as follows. In the remainder of this section we establish 
the notation that will be used throughout the paper. In §2, we present the proper 
orthogonal decomposition and its properties. In §3, we describe the prototypical prob- 
lem used in this article, that is a backward facing step channel flow. We also outline 
the numerical methods used and present numerical results which will later be com- 
pared with reduced order model predictions. In §4, we apply POD for the construction 
of reduced order model. In §5, we formulate an optimal control problem and discuss 
reduced order modeling approach for its solutions. We present computational results 
in §6 with two different control mechanisms. Finally we conclude in §7. 

1.1. Notation. We denote by L 2 (12) the collection of square-int.egrable functions 
defined on flow region 12 C 1R 2 and we denote the associated norm by || • 1 1 o • Let 

H 1 (Q) = Iv £ L 2 (L 2) : — £ L 2 (l 2) for * = 1,2 

and the norm on it be || • 1 1 j . We denote by L 2 (0,T; H l ) the space of all measurable 
functions / : (0, T) — > H 1 such that 


||/||L2(0,T;R>) — 



< oo. 


Vector- valued counterparts of these spaces are denoted by bold-face symbols, e.g., 
= [H 1 (12)] 2 . The L 2 (f 2) or L 2 (12) inner product is denoted by (•, •). We denote 
the inner products for L 2 (r) or L 2 (T) by (•, -)r, where T denotes the boundary of 12. 
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2. The POD Subspace. In order to illustrate the POD subspace and its con- 
struction, we assume for the ease of exposition that we are dealing with the semi- 
discrete nonlinear problem 

^r = £(y,<), ten, y e x, 

where X is a finite dimensional space. If finite element method were used to obtain 
the above semi-discrete problem, X would be a piecewise polynomial space. However, 
the choice for the POD subspace is different. 

2.1. The Proper Orthogonal Decomposition. The underlying problem is 
to identify a structure in a random vector field. Given an ensemble of random vector 
fields U( ! ), we seek to find a function $ which has a structure typical of the members 
of the ensemble. One way to resolve the problem is to project the ensemble on i.e., 
to find $ which is as nearly parallel as possible. Thus we want to maximize 
(U('), $) while removing the amplitude by normalizing it. It is now natural to look 
at a space of functions $ for which the inner-product ($, $) exists, i.e. $ must be 
Z/ 2 (ST) . In order to include the statistics, we must maximize the expression 

in some average sense. Furthermore, since we are only interested in magnitude and 
not the sign, we consider mean of the square of the expression. Following [23], we 
consider ensemble which are “snapshots” , that is the ensemble set 

S = {U« : 1 < i < N } 

are solutions at N different time steps t{ and seek a function $ £ L 2 (f2) that gives 
the best representation of 5 in the sense that it maximizes 

f£LI(u ( E*)| 2 /(*,*). (2-i) 

In other words one seeks a function which has the largest mean square projection on 
the set S. It was shown in [23] that when the number of degrees of freedom required 
to describe U« is larger than the number of snapshots N, it is efficient to express 
the basis function as a linear combination of the snapshots. Thus we assume $ has a 
special form in terms of the original data as 

N 

( 2 . 2 ) 

i = l 

where a; is to determined such that $ maximizes (2.1). The maximization problem 
(2.1) can be cast in an equivalent eigenvalue problem. To see this define, 

i N r 

K <!> = / U (i >(x)U (i >(x')$(x')dx' (2.3) 

i= 1 


(K$,$) = y j U^^(x)$(x)U^^(x / )^(x / )dxdx / 

= ^EK U( ‘U)I 2 . 

*=1 


then 
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Moreover, we have 

(K*,*) _ l( u< 0 .*)l 2 _ . 

(*,*) — (*,*) — A 

Using the calculus of variations, we can find the maximum as described below: Let 
(f>* be a function that maximizes A. We can then write any other function as <j>* +£</>'. 
Then A can be written as 

(K$U$*)+e(K$*,$') + e ( K$/ < $ *)+ e2 ( K$/ - $/ ) _ 

11 ~ ($*,$*) + + e($', $“) + e 2 ($', $') ~ ' 

Clearly, maximum occurs when e = 0 and thus dF Jf’ | e =o = 0. This leads one to 

(K$*,$') = A($*,$')- 

It is now clear that maximization problem (2.1) is the same as finding the eigenvalue 
of the eigenvalue problem 

K$*=A$*. (2.4) 

If (2.2) and (2.3) are introduced into (2.4), we have 

CW = AW, 


where 


( s ) 


Cy= ^/„ u(0(x)uw)(x)rfx 


and W = 


\ w N ) 


It follows from the fact that C is a nonnegative Hermitian matrix that it has a 
complete set of orthogonal eigenvectors 


I !l ) 


( 1 1 ) 


(in 


Wi = 


w 2 = 


W N = 


\ W N ) 


V W N ) 


\<J 


along with a set of eigenvalues Aj > A 2 > > Xn > 0. We can now write down the 

solutions of (2.1): 


N N N 

*!=£>/ U (0 , *2 = X>? u(i) , , = u(i) . 

i — 1 2 = 1 2 = 1 


We also normalize these by requiring 


N 

(w,,w ; ) = y>'-«4* 

2 = 1 


1 

N~\i' 
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It is now easy to check 


(*-«”) = { J ! 


I = m 

l ^ m. 


This completes the construction of the orthonormal set $ 2 , , ^jv}- Then the 

POD subspace is defined as \ POD = span { d> | , d> 2 , , ^jv}. 

To quantify the energy of the data set associated with the corresponding mode 
we note from (2.3) that 




In the next section we will show that POD subspace is optimal in the sense that the 
approximation of the snapshots 




maximizes the captured energy 


£ = iE< u(i) - u(,) ) 


= A* for all Nk < N. 

i = 1 

2.2. Optimality of the basis functions. Given a signal u(x,f) £ L 2 ((0,T) x 
Q) and an approximation u A ' of u in terms of an arbitrary orthonormal basis d', (x) , i =| 

1,2, ,N k : 

N k 

u Ar (x,t) = y 

i — 1 

If ®,-(x) have been nondimensionalized, then the average kinetic energy is given by 


u v ■ u A *c/f2 ) = / cii(t)a* ( t ) \ , 


where (•) denotes the time average operator. We next state a proposition regarding 
the optimality of POD whose proof can be found in [2]. 

Proposition 2 . 1 . Let {$i,$ 2 , denote an orthonormal set of POD 

basis elements and {Aj, , A^v*.} denote the corresponding set of eigenvalues. If 


u A " = 


denotes the approximation to u with respect to the basis, then the following hold: 


(/) ( f 3 i (t)p*(t))=6 1J \ l 


N k N k 


////’) For every N k , (o/?*(o) = Y, Xi > ■ 


i= 1 i = 1 
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where u N (x,t) = J2^=i a;(t)^;(x). 

In essence this proposition states that among all the linear combinations, the one 
corresponds to POD is the best in the sense that it will capture the most kinetic 
energy possible in the average sense. 

2.3. Model reduction aspects. To capture the underlying dynamics of the 
system one needs to keep N sufficiently large. Thus using a Galerkin procedure one 
can obtain a high fidelity model perhaps with large N. However if the eigenvalues of 
the covariance matrix C decays fast, one can choose a cutoff value M <C N and carry 

out a Galerkin procedure with a reduced set of basis elements {d>i, $ 2 , , <&m}- As 

noted earlier, A,- represents the average kinetic energy contained in the first M 

modes. Therefore one can choose M such that A; ~ A through some ex- 

perimentation. Also the ratio — - gives the percentage of the total kinetic energy 

1 * 

contained in the first M POD elements. In fluid flow simulations given below the POD 
system was constructed for N = 100 and the reduced order model was constructed 
with M = 10 which captured 99.99% of the energy. This clearly demonstrates the 
advantage of the reduced order model over the finite element model whose dimension 
was 3,032. 

3. Simulation of Fluid Flow in a Channel. We consider two dimensional 
incompressible fluid flow through a channel with a backward facing step. A schematic 
of the geometry is given in Figure 2. The fluid flow is governed by the Navier-Stokes 
equations which are given by 

u t — -^Au + u ■ Vu + Vp = 0 in Q x (0, T), 

(3.1) 

V ■ u = 0 in ft x (0, T), 

where the velocity u, the pressure p, the time t and the spatial variable x are in 
non-dimensional form. The Reynolds’ number Re is defined as Re = pUoL/p , where 
p is the density, Uq is the nondimensional velocity. The following boundary conditions 
are imposed. 

u = ( 8 ( 2 / — 1/2)(1 — y), 0) = u in on r in x[0,T] 

pn-ikt =( 0 . 0 ) 011 r out x[0,T] 

u =( 0 , 0 ) on r f u r 6 u r s u r c x [o,T], 


The boundary condition on F ov q is not “physical” but used to represent the flow in 
an unbounded region; see [6]. For the finite dimensional approximation and for the 
subsequent reduced order approximation, we need a weak form of the state equations 
(3.1). A weak form of the equations (3.1) has the form; see [27] for similar problems, 


(uf T u Vu, v) + A(Vu, Vv) - (p, V • v) = 0 


(V • u, g) = 0 


(3.2) 


for all test functions (v,g) £ V x Z 2 (ft), where 

V = {veH 1 (fi): v| r \r out = 0}. 

The state variables (u,p) for the problem are taken to be 

u £ Z 2 (0, T; H 1 (f2)), u|r in = u in and u lr\r Qut \r in = 0, 
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p <E L 2 (Q). 

3.1. Finite Dimensional Approximations. To approximate the solutions, we 
will use standard mixed finite element method. For this, we write u and p as linear 
combination of finite number of basis functions: 

N 

— Uq + ^U, :(<)<!> ,:(x), 

2 = 1 

M 

p h = ^2p i (t)'&(x) 

2 = 1 

with Uq be the finite element interpolant of non-zero boundary conditions imposed 
on u. Then 


V A = span{$i, $2 , ffw} 


Q m = spanj'i'i, ® 2 , 

and V jV x Q M C V x L‘ 2 . The approximate system is determined by restricting the 
weak form (3.2) to V' v x Q M with basis functions substituted for the test functions 
v and the basis functions substituted for the test functions q. Then the following 
finite dimensional system results: 


Mu + Su + N(u)u + L 1 p = F, 

Lu = 0, 


(3.3) 


where S is the diffusion matrix, N the convection matrix, L the continuity matrix, 
M the mass matrix and u = Moreover u and p are the finite dimensional 

velocity and pressure, respectively. We call the approximations using standard finite 
element basis functions such as quadratic or linear piecewise polynomials by “full 
order methods/discretization” and the ones using POD by “reduced order methods”. 
For the full discretization we use continuous piecewise quadratics for the velocity 
u and continuous piecewise linear functions for the pressure p; the same triangular 
grid is used for both finite element spaces; This choice of spaces complies with the 
div-stability condition which is required for stable computation of pressure; see [5]. 

The nonlinear differential algebraic equations (DAE) (3.3) is discretized using 
backward Euler in time with the time step At = 0.01 and the resulting nonlinear 
algebraic system is solved using Newtons method along with a banded Gaussian elim- 
ination. 

In Figure 2, the height of the inflow boundary is 0.5 and that of the outflow 
boundary is 1. The length of the narrower section of the channel is 1 and that 
of wider section of the channel is 7. The computational domain was divided into 
682 triangles with finer mesh around the recirculation region. This resulted in a 
system of 3,032 ordinary differential equations that has to be solved for the unknown 
coefficients. We choose throughout in this simulation V max = ^ and H = 1 with a 
corresponding Re = where V),iax = maximum inlet velocity, H =channel height, 
i/=kinematic viscosity of the fluid. 
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It is well known that beyond certain Reynolds’ number the flow separates and a 
recirculation forms near the corner region. We carried out simulations at a Reynolds’ 
number of 1000 and the long term flow simulation is given in Figure 3. It clearly 
predicts the re-circulations first near the corner of the step and the second one near 
the wall opposite to the step. 

4. POD in Flow Simulation. We follow the “snapshots” approach proposed 
in [23] for the derivation of POD basis functions. Let u(x, t ) be a given flow field and 
be the corresponding flow fields at N different time steps tk, he. the 
“snapshots”. We next decompose u(x,f) as follows 

u(x,t) = u m (x) + v(x,t), 

where u m (x) = -h u (x, t*). We also define a spatial correlation matrix C with 

where v ! = v(x,t;). Then the POD basis vectors 4>/,, are defined by 

N 

k = 1 ,N , 

* = 1 

where tu* are the components of the eigenvector W A of the eigenvalue problem 

CW = AW. 

The computation using POD takes the following algorithmic form: 
ALGORITHM I 

(I) Solve the state equation (3.3) at N different time steps and obtain 
“snapshots” S ; See Figure 1. 

(II) Compute the covariant matrix C. The matrix elements of C are 

given by C { j = jr f n v°v j d{2, for i,j = 1,2, , N. 

(III) Solve the eigenvalue problem C W = AW, where C is a nonnegative 

Hermitian matrix and has a complete set of eigenvectors 
Wj.Wj, , Wjy with W l = (u>f,u>5> , w k N ). 

(IV) Obtain the POD basis vectors using <!>,; = Yl!k = i w k wk i 1 < i < N 

and define \ POD = span { 4> i , 4>2, ,4>iv}. And set v = Yli=i a i(t)®i- 

(V) Restrict the weak form (3.2) to \ POD and solve for a,-, i = 1, 2, , N. 

(VI) Set u (x,t) = u m ( x ) +J2? =1 Q '(^) < ^i- 


For the channel flow problem we obtained snapshots u(x,t,) of the flow at 100 
regular intervals. The correlation matrix C was formed with the aid of the finite 
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element routine and the eigenvalue solve was carried out using the RG subroutine in 
the Fortran library EISPACK. The eigenvalue spectrum from the correlation matrix 
C is shown in Figure 4 (left). As shown in the figure, the eigenvalues quickly decay 
and thus very few modes capture the essential energy in the flow. 

4.1. POD Reduced Order Model. In this section, we consider the construc- 
tion of POD reduced order model using a Galerkin projection of the Navier-Stokes 
equations onto a space spanned by the POD basis elements. The nature of the POD 
model is that it requires fewer basis vectors than that is used to approximate the flow 
field. In fact the first M (<C N) modes carry most of the energy in the flow and if we 
choose M such that 


N M 

yi ^ ~ y ’ 

2=1 2=1 

we obtain a reduced order model. We found out 10(= M) POD basis vectors are 
sufficient to capture 99% of the energy. In order to derive the reduced order model, 
let us apply Algorithm I, choose M and expand the solution as 

M 

u(x, t) = u m (x) + ^2 Q 'i(0$i(x). (4-1) 

2 = 1 

Then the Galerkin approximation of the weak form (3.2) is as follows 

(u t + u • Vu, * f ) - (p, V • *j) + V*i) + (pn - = 0, (4.2) 

for all G V POD . At this point it is important to note that the eigenfuntions 
are divergence free as flow is incompressible and satisfy zero boundary conditions on 
r \ F ou ^. Using these properties of and the boundary condition on r ou t, we see 
that the pressure term and the boundary terms vanishes. Then (4.2) reduces to 


(u f + u ■ Vu, * 4 ) + ^(Vu, V*<) = 0, (4.3) 

for all G \ POD . On substitution of (4.1) into (4.3) we obtain the following 
nonlinear evolution equation for the coefficients cq(2): 

q = Aa + a T A r a + e, , 

a(0) = Qo, 1 j 


where 


«<)• — (uo,^), ‘ (iJ-m 

fie 

e»- = -(u m • V u m , $,) - -^-(Vu m ,Vf,j, A fm = ~(®k • V$ ( , $,•). 

He 

The solution to the above initial value problem (4.4) was obtained using an implicit 
Euler method for the coefficients of the POD approximation. 
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4.2. Numerical Results. We selected {<i} at hundred regular time instances 
in the time interval [0, 10]. The eigenvalue analysis of the correlation matrix resulted 
in M=10 and a reduced order system of dimension 10 was constructed. The initial 
value problem (4.4) for the nonlinear ODE was solved using backward Euler method 
with the time step At = 10 -3 and the resulting nonlinear algebraic system was solved 
using Newton iterative method. The Figures 5-6 are the channel flow computations 
with “full solution” and reduced order solution at time t= 10 for various stations in 
the channel which shows excellent qualitative agreement. 

In Table I, we show that £\ norm of the difference between solutions of the POD 
reduced order and full order solution decays as the dimension of the POD subspace 
increases. The percentages of the full order model energy captured by the POD 
reduced order model are also given in Table I which indicates only 9 basis functions 
were enough to capture 99.9% of the energy of the full order model. 

In order to illustrate the features of the POD reduced order model, let us next 
compare it with another reduced order model based on the so called reduced basis 
method (RBM); see [9]. In [9] several ways to choose reduced basis subspaces were 
discussed. Here we consider the so called Lagrange subspace. The basis elements in 
the Lagrange subspace are snapshots of the problem obtained by solving the system 
(3.1) using a full order method. Supposing {'L ! }V. 1 denote the snapshots, the reduced 
order subspace is defined as y RBM = span{'L;}f£ j . Once we have a reduced order 
subspace y RBM , the system (3.1) is projected onto y RBM to obtain a reduced order 
model as in §4.1. 

In algorithmic form the RBM can be summarized in the following form: 

ALGORITHM II 

(I) Solve the state equation (3.3) at N different time steps and obtain 
“snapshots” S\ See Figure 1. 

(II) Set \ RBM = span]^!, <I r 2 , , And set u = u 0 + YiLi a i(t)^it 

where uo account for the nonzero boundary values. 

(III) Restrict the weak form (3.2) to \ RBM and solve for a;, i = 1,2, , M . 

(IV) Set u(x,f) = u 0 (x) + YaLi a i{t)^i- 


For numerical implementation of RBM, we considered, the channel flow case described 
earlier and compared its performance with the full order model. To generate Lagrange 
basis, we obtain snapshots of the model using the full order discretization at M regular 
time instances between t = 0 and t = 10 non-dimensional time. In order to see whether 
the reduced order approximation becomes more accurate as the dimension increases 
we computed the i\ norm of the difference between the reduced order and full order 
solutions. In Table II, we present the t\ norm error using M = 3, 6, 9 and 12 basis 
functions. We also report the condition numbers of the resulting mass matrices. As 
seen, the condition number can increase dramatically with increasing basis elements 
deteriorating convergence. However, the POD reduced order model do not generate 
such bad condition numbers as evidenced in Table I. Moreover, the POD allows easy 
generation of linearly independent basis elements and more stable system matrices. 
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5. All Optimal Control Problem. Minimization of vorticity level in flow do- 
main is of interest in control/delay of transition of flow past bluff bodies. Thus in 
this section we formulate a related optimal control problem in channel flow. Flow 
configuration considered is a backward facing step channel. As Reynolds’ number 
is increased, the flow separates near the corner of the step. The objective of the 
optimal control is to reduce the size of the recirculation and hence of the length of 
re-attachment. The control action is effected either through the movement of a por- 
tion of wall r c or through blowing on T c . In terms of boundary condition it takes the 
following form along the boundary T c , 

u = c(t) g(x) on T c x [0, T], 

where c(t) : [0,T] — > IR and g(x) represent respectively the temporal dependence 
and spatial distribution of the fluid velocity on the boundary T c . The choice of 
cost functional or objective functional to meet the control objective of reducing the 
recirculation is not trivial. Here we will consider a functional of the form 

J{ u) = f 1 1 V x u\\ldt 
Jo 

which corresponds to minimization of vorticity levels in the flow. The task is to find 
c(t) or, rather, its time derivative U = c(t) such that the cost functional 

J(u,U) = - 2 {\\V x u\\l+ i3\U\ 2 } dt 

is minimized subject to the constraints that the flow fields satisfy the Navier-Stokes 
equations. The appearance of the second term in the cost functional J is necessary 
since we will not impose any a priori constraints on the controls. The parameter j3 > 0 
adjusts the relative weight of the two terms in the functional and roughly speaking, 
/? is large for expensive controls and small for inexpensive controls. 

In order to obtain “snapshots” for POD basis functions, we introduce 

u c (.e) = u Cl (x) - u Co (x), 

where u Cl is a steady flow with c(i, ) = 0.1 on r c and u c ° is that with c(t) = 0 on F c . 
Then the “snapshots” are defined as 

u{x,t k ) - c(t k ) u c (x) 

and the basis functions <5; as defined in Algorithm I have zero boundary conditions 
on the Dirichlet boundaries. The velocity expansion is defined as 

M 

u(x, /) = u m (x) + c(f)u c (x) + ^2 oti(t)$i(x) (5.1) 

2 = 1 

so as to automatically satisfies all the Dirichlet boundary conditions. 

5.1. The Reduced Order Control Problem. Inserting the expansion (5.1) 
into the Galerkin projection (4.3) of the Navier-Stokes equations, we obtain 

Otf T An T €t Not T U a T (b -}- bci)c -}- c d -(- e = 0, 


(5.2) 
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where 


ii = (u c , <!>,■), di = (u c ■ Vu c , $;), 


bi = (u m • Vu c , $,■) + (u c • Vu m , $;) + — (Vu c ,V$i), 

tie 


Bij = ($/ • Vu c , $;) + (u c • V$j. $,;) 


and .4, A/" and e are as defined in §4. Setting 



^ a T Not + (b + Ba)c + c 2 d + <3 ^ 


we obtain the reduced order control problem: 

rT 


Minimize J(X,U) 


n 


B 


L(X) + iu 2 \ dt 


subject to 


X = F(X) + BU, 
X(0) = X o , 


(5.3) 


(5.4) 


where 

L(X) = ]-X T QX + X-f 1 + f 2 , and F(X) = -AX - N(X). 

At this point one can employ a variety of numerical methods designed for finite dimen- 
sional nonlinear optimal control problems such as multiple shooting methods. Our 
method here is based on Newtons method for the necessary condition of optimality 
or the so called optimality system albeit for the discrete version of it. 

We further remark here that finite dimensional control systems like the one given 
above can also be obtained using, for example, finite element method as in §3. How- 
ever, their size is too large for practical control systems whereas POD based reduced 
order control systems are low order and maintain high fidelity. This makes our ap- 
proach extremely attractive for optimal control problems governed by partial differ- 
ential equations. 

5.2. Approximation of the Reduced Order Control Problem. We con- 
sider the second order discrete time approximation of (5.3)-(5.4) obtained by using 
the Crank-Nicholson for time discretization of (5.4) and the trapezoidal rule for the 
discretization of the integral in the cost functional (5.3). We obtain 


Minimize J K = eU (L(X k ~ 1 ) + L(X k )) + -\U k \A At 


(5.5) 
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subject to 

^—^— = ^(F(X k ) + F{X k ~ 1 ))+BU k , 1 <k<K (5.6) 

where K At = T. Note that if uiAt < 1 then the mapping <f>(X) = X — A tF(X) is 
dissipative, that is 

(F(X!) - F(X 2 ),Xi - X 2 ) > (1 - uAt) |Xi - X 2 | 2 . 

Thus for U = {U k }J? =1 , there exists a unique X = {X fc } A =1 satisfying the constraint 
(5.6) and depending continuously on U . Moreover if wAt < 1 then there exists an 
optimal pair (X k , U k ) to problem (5.5)-(5.6). 

The necessary optimality condition, see Appendix A, for (5.5)-(5.6) is given by 

' X " A f | (^(X fc ) + F(X»-*)) - 1 BB V 

(5.7) 

_ - At C = 2 F X (*. k f (C k +C k+1 ) + L x (X k ) 

for 1 < k < K, with X° = Xo and £ A+1 = 0, and the optimal control U k to (5.5)— (5.6) 
is given by 

U k = ~^ B T C fc , \<k< K. (5.8) 

Assume that {U m }^ =1 is a sequence of solutions to (5.5) — (5.6) with associated state 
and adjoint states {(X m ,£ m )}~ =1 such that (5.7) holds. Let U m denote the step 
function defined by U m (t) = U k on (tk-\ , t k ), 1 < k < K and X m and £ be the 
piecewise linear functions defined by 


X m (t) =X te ~ 1 + 


-X 

"a T 


fe-i 


(t-tk- 1 ) and C m (t ) = ( k + ~ 


fc + 1 


At 


{t - t k _ i). 


~ ~ ~m 

Then it can be proved that the sequence (X m , U m X ) has a convergent subsequence 
as At — ► 0 and for every cluster point (X, U, Q, U is an optimal control of (5.3)-(5.4) 
and (X,U,Q satisfies the necessary optimality condition: 


f X(t) = F(X(t)) - 
l -C = F x (X) T C(t) + L x (X) 


(5.9) 


with X(0) = Xo and £(T) = 0. Furthermore, (5.7) is a sparse system of nonlinear 
equations for (X, £). We solve the system (5.7) by Newton method in our calculations 
of the optimal control {U k }. The Jacobian J of equation (5.7) has the sparse structure; 


J = 


A 

Q 


s 

-A f 


where S and Q are block-wise diagonal and A is block-wise lower bi-diagonal with 
block size K. The diagonal block A kk and the off-diagonal block A^ + i k of A are 
given by 

Ak ’ k = jAt + \ Fx(X " } and Ak+1 ’ k = “2 L + l Fx(Xk) - 
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S has the constant diagonal block ^ B r B and the diagonal block Qk k of Q is given 

by 


Qk, k = L xx (X. k ) + J2 (^)x.x(X fc ). 

i = 1 

6. Numerical Results. Here we present numerical results for the POD based 
control and compare its performance with that of RBM. The flow configuration is 
chosen as the two dimensional backward facing step. The control objective is to 
reduce the recirculation behind the step and thus of the re-attachment length. The 
cost functional is taken to be the vorticity functional defined earlier. Two control 
mechanisms are considered. In the first example, we consider moving wall as control 
and in the second we consider blowing of mass through a portion of the boundary. 

6.1. Example I. In this example, the control is introduced into the problem 
through the movement of wall on the boundary T c : 

u = c(t)r on r c , 
where r is the unit tangential vector. 

6.1.1. Test I (POD). We present numerical results for POD approach in solv- 
ing the optimal control problem at Re = 200. Recall that the control problem we 
consider is 

Minimize J(u, U ) = \ fj { || V x u||q + /?|f/| 2 } dt 
subject to 


(uf T u • Vu, T,;) + ^(Vu, Vd>i) = 0, i=l,....,M, 


where M is the number of POD modes and 

M 

u = u m + c(t) u c + y^o-;(t)<f>,-(x). 

2 = 1 

The initial conditions for the states and controls are u (x, 0) = 0 and c(0) = 0, 
respectively. Our choice of the portion of the boundary r c , where control is applied, 
is the line segment between y = 0 and y = 0.5 at x = 1; see Figure 2. This choice 
here is motivated by the fact that if one wants maximum influence in the flow, then 
the control has to be in that vicinity. The time interval was chosen [0, T] with T = 10 
and the number of POD modes was taken as 4. In Figures 7-9 we present numerical 
results for the penalty parameter /? = -gb and the time step At = 0.01. The numerical 
solution of the optimal control problem was obtained using the Newton’s method 
described in §5.2. The flow fields presented in Figures 7-8 are u component of the 
flow field u for the controlled and baseline cases at different stations in the channel. 
As one can see when control is applied the u velocity becomes positive where it is 
otherwise negative and helps the formation of recirculation. Significant reduction in 
the recirculation bubble and re-attachment length were also observed compared to the 
uncontrolled case. We also carried out our calculation with different initial conditions 
Xq and the results were qualitatively similar to those described above. 
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6.1.2. Test II (RBM). We present here numerical results for the RBM ap- 
proach to the same optimal control problem as in Test I. Like in POD, the solution 
expansion in RBM is of the form 


M M 0 

11 = lip + ai(t)<Pi + '^2 c s -(t)ui 
i = 1 i—1 

where uo denotes the flow corresponding to a zero control, i.e, it satisfies uo = 0 on 

r c , and ui, ,ujf 0 denotes solution of (3.1) with nonzero values on the control 

part of the boundary T c . For the numerical results presented in Figure 9, we have 
used basis functions based on the data shown in Table III. For the simulation of the 
control problem we take 


M 

u = u 0 +^ «;(<)$; + c{t)$ M + 1, 

2 — 1 

where M = 4 and the basis functions d>i , d>2, $3, d>4 and <J>5 are chosen as <I>i = 
u 4 — 2ui + uo, $2 = U3 — 3ui + 2u 0 , <I>3 = u 4 — 4ii! + 3u 0 , <J> 4 = u 5 — 5ui + 4uo and 
$ 5 = uo — ui. The time interval [0,T], the time step At, the penalty parameter and 
the other date were all taken the same as in the previous test case. The numerical 
solution of the control problem was also computed using the same method. The control 
distribution presented in Figure 9 and the controlled flow fields (not presented here) 
all agree well with that of POD. This shows the ability of RBM to provide very good 
controls with very few elements. However, RBM can be sensitive in terms of condition 
numbers of the system matrices as one increases the number of basis functions in order 
to improve convergence and accuracy. 

6.2. Example II. In this example, we consider a different control mechanism 
from the previous one. The control is effected through blowing on the lower quarter 
of the boundary T c . Thus we consider 

f c(2)g(x) on 0 <y<\ 

l 0 on \<y <\ 

and g(x) = (10j/(g — y),0). The initial conditions for the states and control were 
u(x, 0) = 0 and c(0) = 0. 

6.2.1. Test I (POD). We present numerical results for POD approach in solv- 
ing the optimal control problem at Re = 500. The control is introduced into the 
problem through the blowing on the lower quarter of the boundary T c . The compu- 
tational domain was similar to Figure 2, but the length of the narrower section of 
the channel is 0.5 and that of wider section of the channel is 12. The computational 
domain was divided into 794 triangles with finer mesh around the recirculation re- 
gion. Our choice of the portion of the boundary, where control is applied, is the line 
segment between y — 0 and y = 0.125 at x = 0.5. A linear time varying profile for c(t) 
was used to generate 100 “snapshots” for the generation of POD modes. We carried 
out calculations with 4, 9 and 14 modes in the time interval [0,T] with T — 10. In 
Figures 10-13, we present numerical results for the penalty parameter /? = gh The 
numerical solution of the optimal control problem was obtained using the Newton’s 
method described in §5.2. The initial conditions for the state and adjoint state were 
all zero and the convergence tolerance was taken to be 10“ '. 



16 


S.S. RAVINDRAN 


The controlled flow fields with 4, 9 and 14 modes showed similar results and 
hence only results with 9 modes are presented. The flow fields presented in Figures 
10-11 are u component of the flow field u at different stations in the channel for the 
controlled and uncontrolled cases with 9 POD modes. As indicated by the controlled 
flow fields, separation has been effectively eliminated by the optimal blowing control. 
Significant reduction in the recirculation bubble and re-attachment length were also 
observed compared to the uncontrolled case. 

6.2.2. Test II (RBM). We present here numerical results for the RBM ap- 
proach. Like in POD, the solution expansion in RBM is of the form 

M M 0 

U = U 0 + ^2 + ^2 c *'(*) u * 

2=1 2—1 

where uo denotes the flow corresponding to a zero control, i.e, it satisfies uo = 0 on 

r c , and ui, j U M 0 denotes solution of (3.1) with nonzero values on the control 

part of the boundary T c . For the numerical results presented in Figure 13, we have 
used basis functions based on the data shown in Table IV. For the simulation of the 
control problem we take 

M 

u = $i + (*)$,• + c(t)$ M + 1, 

2 — 1 

with M=4, 9 and 14, respectively. The RBM basis function selection was simi- 
lar to that of the previous example. For example, when M=4, the basis functions 
,$2j $3, d>4. $5 and <E>6 were chosen as <f>i = ui — 0 . 1 (rig — ui)/9.9, $2 = 112 — 
1.9(u 6 — ui)/9.9 — ui, $3 = u 3 — 3.9(u 6 — ui)/9.9 — iii, $4 = ii5-5.9(u 6 -ui)/9.9-ui, 
= U6 — 7.9 (u 6 — ui)/9.9 — ui and $g = ug — Ui . The time interval [0, T], the time 
step At, the penalty parameter and the other data were all taken the same as in the 
previous case. The numerical solution of the control problem was also computed using 
the same method. The control distribution presented in Figure 13 and the controlled 
flow fields (not presented here) all agree well with that of POD. These results seem 
to re-confirm our earlier observations about RBM’s performance and the effectiveness 
of the control mechanism. 

7. Conclusion. In this article we have presented a reduced order modeling ap- 
proach for optimal control of fluid flows. The reduced order models suitable for control 
and which captures the essential physics were developed using the POD. Our com- 
putational investigations into the use of reduced order methods for control suggest 
promise. Significant computational savings were evidenced in the test cases consid- 
ered. In the the reduced basis method there is no systematic way to increase the 
level of accuracy, and ill-conditioned system matrices can make it impossible to im- 
prove the accuracy. However, the POD provides a systematic and optimal way to 
improve the level of accuracy while maintaining well-conditioned system matrices. As 
we can see, these are not provided as a generic methods rather they must be used 
with care. Whenever they can be effective they can provide significant performance 
with substantially lower on-line computational resources. We have also investigated 
the feasibility of two different control mechanisms. Controlled flow fields in both cases 
were comparable. This seems to indicate that the choice of the type of control (wall 
movement or blowing) should be decided based on which type is the easiest to apply 
in the particular application. 
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Appendix A. We consider the problem of minimizing J(U) in (5.3) subject to 
the ordinary differential equations and initial conditions in (5.4). The control U* is 
extremal and J has a relative minimum, if there exists an e such that for all functions 
satisfying \\U — U*\\ < e the difference J(U) — J(U*) > 0. We have the following 
classical theorem. 

Theorem A.l. For U* to be extremal, it is necessary that SJ(U*,SU*) = 0 for 
all SU , where SJ is the variation in J with respect to the variation in SU in U . 

A proof of this theorem can be found in [13], In order to apply this theorem, let 
us introduce a vector of Lagrange multipliers 

C T = (Ci, Xn) 

and form an augmented functional including the constraints 

J = f Q (i(X) + (X - /(X) - B U) dt. (.4.1) 

We integrate by parts and take variations in J corresponding to variations SU in U 
to get 

SJ = [-C T SX] t=T + J |ix(X) + I3USU+C SX + C t F x 6X + C t BSu\ dt, (.4.2) 

note that 4X(0) = 0 as X(0) is given. We may eliminate some of the terms in (A. 2) 
by defining 

C=-Lx-C T F(X), and C(0) = 0. (A3) 

Equation (A. 3) then reduces to 

pT 

SJ = f [/3U+£B}SUdt 
J o 

and now from the Theorem A.l, a necessary condition for U* to be extremal is that 

fdU + C T B = 0. 

The state equation (5.4) and the adjoint equation (A. 3) form 2 n differential equations 
with boundary conditions X(0) = Xo and £{T) = 0. 

We like to remark here that the Theorem A.l provides a necessary condition for 
optimality. It is easy to show the fact that it is in general not a sufficient optimality 
condition, i.e. (X* ,U* X*) can be an extremal element without (X*,U*) being a 
solution to (5.3)— (5.4) . 
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TABLE I. l\ norm difference between full order and POD reduced order model 
solutions, condition numbers of the mass matrix and percentage of full order model 
energy captured with M = 3,6, 9, 12, 15 and 20 POD basis functions. 


M 

3 

6 

9 

12 

15 

20 

l\ Error 

0.0013 

0.001 

.000769 

0.000359 

0.00029 

0.00017 

Condition k 

1.0 

1.0 

1.0 

1.0 

1.0 

1.0 

% of Energy 

97.0 

99.68 

99.96 

99.997 

99.999 

99.9999 


TABLE II. l\ norm difference between full order and. RBM reduced order model 
solutions and condition numbers of the mass matrix with M = 3,6,9 and 12 basis 
functions. 


M 

3 

6 

9 

12 

I\ Error 

0.0327 

0.0057 

.0035 

0.0023 

Condition k 

613.45 

21840 

388556 

1974595 


TABLE III. Simulation data for RBM reduced order model; Example I. 



u 0 

Ui 

u 2 

u 3 

11/} 

u 5 

Control 

0 

-0.1 

-0.2 

-0.3 

-0.4 

-0.5 


TABLE IV. Simulation data for RBM reduced order model with f, 9 and If mode, 
respectively; Example II. 



Ul 

u 2 

u 3 

u 4 

Us 

u 6 

U7 

U 8 

U 9 

UlO 

Ull 

Ul2 

Ul3 

Ul4 

Ul5 

Uie 1 

Control 

0.1 

2 

4 

mm 

8 

El 











Control 

0.1 

1 

2 

mm 

4 

mm 

6 

7 

8 

9 

10 





1 

Control 

0.1 

1 

2 

3 

4 

5 

5.5 

6 

6.5 

7 

7.5 

8 

8.5 

9 

9.5 

10 l 
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FIG. 1. Illustration of sample selection (“snapshots”) for reduced order methods. 



FIG. 2. Computational domain for the backward-facing-step channel problem 



FIG. 3. Velocity field at t = 10 and Re=1000. 
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FIG. 4. Eigenvalues of the correlation matrix for the baseline and controlled cases 
(Test I and II), respectively 



FIG. 5. Comparison of full model versus reduced model solutions; velocity compo- 
nents u and v at x = 1.31, x = 1.9 and x = 2.38 



FIG. 6. Comparison of full model versus reduced model solutions; velocity com- 
ponents u and v at x = 2.53, x = 2.69 and x = 2.84 
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